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ABSTRACT 

Observed galaxy clusters often exhibit X-ray morphologies suggestive of recent interaction 
with an infalling subcluster. Abell 3376 is a nearby (z — 0.046) massive galaxy cluster whose 
bullet-shaped X-ray emission indicates that it may have undergone a recent collision. It dis- 
plays a pair of Mpc-scale radio relics and its brightest cluster galaxy is located 970hy^ kpc 
away from the peak of X-ray emission, where the second brightest galaxy lies. We attempt 
to recover the dynamical history of Abell 3376. We perform a set of A-body adiabatic hy- 
drodynamical simulations using the SPH code Gadget-2. These simulations of binary cluster 
collisions are aimed at exploring the parameter space of possible initial configurations. By 
attempting to match X-ray morphology, temperature, virial mass and X-ray luminosity, we 
set approximate constraints on some merger parameters. Our best models suggest a collision 
of clusters with mass ratio in the range 1/6-1/8, and having a subcluster with central gas 
density four times higher than that of the major cluster. Models with small impact parameter 
(6 < 150 kpc), if any, are preferred. We estimate that Abell 3376 is observed approximately 
0.5 Gyr after core passage, and that the collision axis is inclined by i « 40° with respect to 
the plane of the sky. The infalling subcluster drives a supersonic shock wave that propagates 
at almost 2600 km/s, implying a Mach number as high as M ~ 4; but we show how it would 
have been underestimated as M. ~ 3 due to projection effects. 
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1 INTRODUCTION 

Mergers are the processes through which galaxy clusters assemble 
in the hierarchical scenario of structure formation. Observed galaxy 
clusters often display perturbed morphologies suggesting they un- 
derwent recent interactions with less massive subclusters. 

A number of galaxy clusters exhibit diffuse radio emission in 
their periphery, which is not associated with galactic sources (e.g. 



Feretti & Giovanninill99f3:lvan Weeren et al.l2009l;lBonafede et al.l 
2012b . These Mpc-scale structures are know as radio relics and they 
are useful as probes of merger shocks. They are generally inter- 
preted as being driven by shock waves that propagate outwards at 
supersonic speeds, re-accelerating relativistic electrons through the 
Fer mi mechanism in the very low-density o utskirts of the cluster 
(e.g. lFuiita et alj2003l ; lGa"bici & Blasil2003h . 

From the theoretical standpoint, idealised numerical sim- 
ulations of binary cluster mergers suppl y a wealth of i nsight 
into the outcomes of these events (e.g . Roettiger et al.l 1 19931: 



Schindler & Mueller 1993 ; IPearce et al.l Il994j; IRoettiger et alj 
19971 : iRickerl Il998l: IRoettiger & Floresl 120001: iRicker & SarazirJ 
200l |; iRitchie & Thomas! |2002| ; IPoole et alj l200d : IZuHone et at] 
20091 |2010| : IZuHond 1201 If) . Hydrodynamical simulations of 
merging cluster designed to model specific observed objects 

* E-mail: rgmachado@astro.iag.usp.br 



have freque n tly focused on the Bulle t Cluster (e.g. iTakizawal 
20051 120061: iMilosavlievic et alj 120071 : [Springel & Farrarl |2007| : 
Mastropietro & Burkertl l2008l) . Recently. Ivan Weeren et al.l d201 lb 
carried out hydrodynamical simulations to model the galaxy clus- 
ter CIZA J2242. 8+5301 and used its double radio relics to set 
constraints on the m erger geometry, mass ratio and time-scale. 
iBriiggen et aT] d2012b modelled 1RXS J0603. 3+4214, another ra- 
dio relic cluster, as a triple merger. 

In fully cosmological hydrodynamical simulations, clusters 
mergers are studied in a more realistic but less well-controlled 
environment, and sometimes at the cost of lower spatial or mass 
resolution. They provide useful analyses of the statistical proper- 
ties of merger shocks, such as cold fronts (Hallman et al. [2010i) 
and the distributions o f Mach numbers (e.g. | Vazza et al. 1201 lL 
lArava-Melo et al.1 120121 : IPlanelles & Ouilisl I2012T) . Typical Mach 
numbers in shocks driven by cluster collisions tend to be < 3, 
but strong shocks may arise under some circumstances. For exam- 
ple, iFinoguenov et al.1 d2010h obtained A 4 ~ 2 from the density 
and temperature drops in Abell 3667 and iMarkevitch et"ai1 ( 2005 ) 
had found a similar valu e for Abell 520. A sho ck as strong as 
M ~ 4 was estimated by Ivan Weeren et at] d2010b from the radio 
relics of CIZA J2242.8+5301 . 

Using Suzaku X-ray observations, Akamatsu et al. I d2012b 
measured the temperature jump in the western radio relic of 
Abell 3376, obtaining M = 2.91 ± 0.91. From polarisation 
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and spectral indices s tudies of the radio relics of Abell 3376, 
iKale et alj l201ll . l2012h estimated M ~ 2.2 - 3.3. From a cos- 
mological simulation. [Paul et alj d201ll) were able to identify one 
merging cluster whose shock structure resembles the radio relics of 
Abell 3376. 

Abell 3376 (hereafter A3376) is a nearby (z = 0.046) massive 
galaxy cluster with a distinctive bullet-like morphology, suggesting 
an ongoing collision taking place. It is possibly the closest cluster 
exhibiting such a morphology. In its outskirts, a pair o f prominent 
Mpc-scale radio relics are present (Bagchi et al. 2006|). Its virial 
mass is estimated as 5.0 x 10 14 Mm dGirardi et al.lll998l based on 
galaxy velocity dis persion) and its 0. 1- 2.4 keV luminosity is Lx — 
2.5 x 10 44 erg/s dEbeling et all 1 19961 based on Rosat data). The 
first and second brightest cluster galaxies (BCG) are separated by 
~ 970/i^q 1 kpc, but it is the second BCG that coincides with the 
pea k of X-ray emission. 

iBagchi et al.1 12006J) have analysed early XMM-Newton X-ray 
and VLA radio data, proposing that the radio relics observed may 
also be the site of acceleration of very energetic cosmic rays, up to 
10 18 eV. This would be possible with first-order Fermi mechanism 
at the Shockwave front. They set forth two possible scenarios that 
could be responsible for producing the shock: (i) the accretion of 
intergalactic medium flowing down towards the cluster; or (ii) the 
collision of two clusters, but the issue remained unresolved. 

Our goal in this work is to run A^-body+SPH simulations 
tailored for the specific case of A3376 using more recent XMM- 
Newton data in order to constrain the actual dynamical history of 
this system. 

This paper is organized as follows. Simulation techniques and 
initial conditions are described in Section|2] X-ray observations are 
described in Section [3] In Section [4] we present the global merger 
evolution, explore the parameter space of different possible col- 
lisions, and compare the simulation results to observations; Mach 
number and dark matter distribution are discussed. Finally, we sum- 
marise and conclude in Section|5] Throughout this work we assume 
a standard ACDM cosmology with Aa = 0.7, Q,m = 0.3 and 
Jfo = 70kms _1 Mpc -1 . 



2 SIMULATIONS 

We set up an idealised numerical model to represent the merg- 
ing of two initially isolated galaxy clusters. The main puipose of 
these simulations is to reproduce certain features of the galaxy 
cluster A3376, especially the global morphology of the intracluster 
medium. Even though this model relies on several simplifications, 
it allows us to reconstruct a possible scenario for the dynamical his- 
tory of A3376. By comparing the results of a large set of simula- 
tions, it is possible to set approximate constraints on some collision 
parameters. 



2.1 Techniques 

We consider the collision of two spherically symmetric galaxy clus- 
ters: a more massive cluster A (the major cluster) and a less mas- 
sive cluster B (also referred to as the minor cluster or the subclus- 
ter). A range of mass ratios Mb /Ma is explored. Each cluster is 
composed of dark matter particles and gas particles, and a baryon 
fraction of 0.18 is adopted throughout. 

In all simulations, the major cluster has a total of 
Na = 6 x 10 6 particles proportionally divided between dark 
matter and gas, such that the mass resolution is the same 



m, = 1 x 10 Mq for both species. The minor cluster's mass and 
particle number are scaled down such that all particles have the 
same mass. 

Simulations were perf ormed using th e public version of the 
parallel SPH code Gadget-2 JSpringefeOoH) with a softening length 
of e = 5 kpc and a maximum time step of 0.001 Gyr. The intra- 
cluster medium is represented by an ideal gas with adiabatic index 
7 = 5/3. As a first approximation, cooling is not taken into account 
in the simulations, since the cooling time-scale is larger than the 
merger time-scale. Galaxies themselves a ccount for a smal l frac- 
tion of the total cluster mass (~ 3%, e.g. lLagana et~ai] [2008) and 
their gravitational influence may be disregarded. Since we do not 
include stars in these simulations, feedback and star formation are 
not present. Even though magnetic fields are believed to give rise to 
the radio relics observed in A3376, they are not expected to play an 
important role in determining the global cluster morphology, and 
are ignored in our simulations. The evolution of the system is fol- 
lowed for 5 Gyr, but the relevant phases take place in a time-scale 
of not more than about 1 Gyr. Because of that, and of the small 
spatial extent of the system, cosmological expansion is neglected. 
Simulations were carried out on a 2304-core SGI Altix cluster. 



2.2 Initial conditions and models 

For the dark matter halo, a iHernauistl dl990h density profile is 
adopted: 
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where Mh is the dark matter halo mass, and ru is a scale length . 
This is similar to the NFW profile (Navarro. Frenk & Whitell 19971) 
except in the outermost parts beyond the r2oo radius (also under- 
stood as the virial radius), with the advantage of having a finite 
total mass. Conveniently, many of its properties (such as potential, 
cumulative mass and distribution function) can be expressed ana- 
lytically. 

For the gas distribution, we employ a iDehnenl Jl993h density 
profile: 
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where M g is the gas mass and r g is a scale length. This has the ben- 
efit of preserving the analytical simplicity of several useful quan- 
tities (chiefly the derivatives of the potential), while also allowing 
the possibility of a flat core, which is achieved by setting its 7 pa- 
rameter to zero. The resulting profi le resembles that of a /3-model 
dCavaliere & Fusco-Femian ol 19761) and is thus suitable to represent 
the intracluster medium of an undisturbed cluster without a pro- 
nounced cool-core, i.e. without a steep density profile in the centre. 

To create a numerical realisation of this dark matter + gas sys- 
tem, we set up the i nitial positions and veloci ties following the pro- 
cedure outlined by Kazant zidis et alj j2006l) . First, the dark mat- 
ter cumulative mass function is uniformly sampled in the interval 
[0, Mh] and the inverse function r(M) is used to provide r, the 
distance from the centre, for each dark matter particle. The same is 
done for the gas particles. From these radii, cartesian coordinates 
are obtained by assigning random directions to the position vectors. 

One possible method to obtain the velocities of the colli- 
sionless particles is the so called 'local Maxwellian approxima- 
tion'. It consists in drawing velocities from Maxwellian distribu- 
tions with dispe rsions obtained from solv ing the Jeans equation 
at each radius dBinnev & Tremainel fl987). However, due to the 
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inadequacies of this a pproach which have been pointed out by 
iKazantzidis et alj J2004 ). one must resort to the distribution func- 
tion itself. In this way, no assumptions need to be made about the 
local shape of the velocity distribution (apart from the assump- 
tions of isotropy and spherical symmetry) and the exact distribution 
function f(£) is given by Eddington's formula ( lEddingtodll916l: 
iBinnev & Tremaine|[l987l) : 



Table 1. Inital condition parameters. Model 233 is the fiducial model. 
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where 4< = — $ = — + <I>g) is the relative total potential, 
and £ — \I/ — u 2 /2 is the relative energy. For physically mean- 
ingful Ph(r) and ^(r), the last term in equation ((3} vanishes. This 
leaves an integrand that depends on (the second derivative of) the 
dark matter density expressed as a function of the total potential, 
ph = ph(^)- The integration is evaluated numerically and the 
function /(£) is tabulated in a fine grid over a range of energies 
and then interpolated wherever necessary. Random pairs (£, f) are 
drawn and values of v 2 are accepted according to the I von Neumann! 
( 1951) rejection technique. This assigns a speed to each collision- 
less particle and cartesian velocity components are obtained assum- 
ing random directions for the velocity vectors. 

Since the gas is assumed to be in hydrostatic equilibrium, each 
volume element of the fluid is initially at rest. The SPH particles 
require an additional quantity to be set up: their internal energy (i.e. 
temperature). From the assumption of hydrostatic equilibrium, it 
follows that for a chosen gas density p g (r), the temperature profile 
is uniquely specified as: 
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where fi is the mean molecular weight, m p is the proton mass, k 
is the Boltzmann constant and M(r') = M g (r') + M h (r') is the 
total mass inside the radius r'. 

Each cluster model is allowed to relax in isolation for a period 
of 5 Gyr, typically a few dynamical time-scales, prior to the begin- 
ning of the actual collision. This ensures that transient numerical 
effects, however minor, will have had time to subside. The clusters 
are then placed at a sufficiently large initial separation do, having 
an initial relative velocity Vo in the direction of the x-axis. To avoid 
spurious tidal effects in the initial conditions, we employ an initial 
separation of do = 4 Mpc, which is approximately twice as large 
as the sum of the clusters' virial radii. 

Numerous combinations of plausible initial condition param- 
eters are possible and, in the search for a 'best-fitting' model, hun- 
dreds of simulations were run. Table [T] lists the initial condition 
parameters of the sample of models that are reported in this paper. 
This sample focuses on the variations of four parameters, namely: 
the total mass ratio Mb /Ma, ranging from 1/2 to 1/8; the ratio of 
the central gas densities no,s/«-o,A ranging from 1 to 6; the initial 
relative velocity vo, ranging from 500 to 2000 km/s; and the impact 
parameter bo, ranging from to 500 kpc. 

The major cluster's total mass is always M A = 6 x 1O 14 M 
and its central gas density is no, a = 1.2 x 10~ 2 cm~ 3 in all cases. 
These values were chosen to ensure that both the total mass and the 
total X-ray luminosity of the resulting object are within the same 
order of magnitude as the observed cluster. The initial relative ve- 
locity is always parallel to the s-axis, even when a non-zero impact 
parameter (a shift in the initial position of the subcluster along the 
y-axis direction) is present. 
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3 X-RAY DATA 
3.1 Observations 

The numerical simulations presented here will be compared to 
archival X-ray observations. A3376 was observed twice by the 
XMM-Newton satellite, in 2003 (revolution 0606, RI. M. Marke- 
vitch) and in 2007 (revolution 1411, RI. M. Johns ton-Hollitt). 
While the 2003 observation was partially analysed by Bagc hi et al.1 
(2006) (they did not use the pn detector), the 2007 observation is, 
as far as we know, unpublished. 

Both observations were done in Prime Full Window with 
"medium" filter. We have run the Science Analysis System (SASQ 
11.0) pipeline, removing bad pixels, electronic noise, and correct- 
ing for charge transfer losses. For the EPIC MOS 1 and MOS2 cam- 
eras we have kept only events with PATTERN 12 and FLAG = 
(events on the field of view). For the pn camera, have kept events 
with PATTERN 4 and FLAG = 0, following the standard pro- 
cedure recommended by the SAS team. Both observations were 
screened for high particle background periods. We have constructed 
light-curves in the [1-12 keV] band and filtered out time intervals 
of anomalously high flux. The final exposure times for the 2003 ob- 
servation were 23.0, 22.9, and 16.0 ks for the MOS1, MOS2, and 
pn, respectively. For the 2007 observation they were 36.8, 39.8, and 
25.6 ks for the MOS1, MOS2, and pn, respectively. 

When necessary, the background was taken into account by 
using the publicly avail able EPIC blank sky templates described by 
iRead & Ponmanl 120031) . taking into account the observation mode 
and filter. Each blank sky background was further normalised using 
the observed spectrum obtained in an annulus (between 12.5-14.0 
arcmin), taking care to avoid the cluster emission. 



3.2 Analysis 

With the cleaned event files, we have produced the exposure-map 
corrected images in the [0.5-8.0 keV] band combining all XMM 
data available. This image will be compared below to a simulated 
image from the TV-body simulation. 



1 See http://xmm.esac.esa.int/ 
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We have also pro duced a 2D temperature map , using an adap- 
tive kernel technique fourret & Lima Netoll2008l) . A cell of vari- 
able size must have a minimum count number (of the order of 
10 3 after background subtraction). A cell that meets this criterium 
has its spectrum fitted by an absorbed, singl e temperature MEKAL 
model (bremsstrahlu ng and emission lines. iKaastra & Mewdll993l; 
Liedahletal]|l995h using the XSPEC 12.5 packagj 2 ! lArnaudl 
19961) . The free parameters are the intra-cluster plasma tempera- 



ture and metallicity (metal abundance, mainly iron). 

The spectral fits are done in the [0.7-8.0 keV] band. We fixed 
the absorption using the Galactic value of the neutral hydrogen col- 
umn density (4.6 x 10 20 cm -2 ; Leiden/ Argentine/Bonn (LAB) Sur- 
vey), estimated with the nh task from FTOOLSjf). The effective area 
files (ARFs) and the response matrices (RMFs) were computed for 
each cell in the image grid. We produced a temperature map com- 
bining both observations, which will be compared to our simula- 
tions below. 



4 RESULTS 

4.1 Global evolution 

To illustrate the temporal evolution of a cluster merger simulation, 
Fig. [U displays a sequence of snapshots (of model 233) spaced by 
0.2 Gyr, in frames of 1.6x1.6 Mpc. This is a head-on collision 
along the a;-axis, which appears arbitrarily rotated on the plane of 
the sky merely to match the position angle of the observation, a 
rotation having of course no further relevance. The meaningful ori- 
entation is the inclination angle i, between the collision axis and the 
plane of the sky. All frames in Fig.Q]are projected at a inclination 
i = 40°, discussed in detail in Section l4.2.5l To allow for a some- 
what fair comparison with observations, we compute the projected 
X-ray surface brightness (left column of Fig.Q] and the projected 
emission- weighted temperature (right column of Fig. [T}. 

Even though the simulations themselves are adiabatic, when 
generating these simulated images we assume the X-ray radia- 
tive losses can be described by a cooling function that takes more 
than just bremsstrahlung into account. From the simulation output, 
we estimate the emission using a cooling function A(T) that was 
computed using the MEKAL emission spectrum model with the 
XSPEC 12.5 package, appropriate for a plasma of metallicity equal 
to 0.3 Zq. At high energies the free-free emission (oc n 2 \/T) is 
predominant (and that is the quantity usually employed as a proxy 
for emission in galaxy cluster simulations), but at low energies col- 
lisional excitation dominates. The emission is then projected along 
the line of sight to give the X-ray surface brightness maps. Tem- 
perature maps are generated by weighting particle temperatures by 
their emission, and projecting them along the line of sight. 

In Fig. [T] the subcluster comes from the lower right corner of 
the frame towards the major cluster. Most of the relevant dynamical 
evolution of the merger takes place within an interval of 1 Gyr. The 
simulation starts at t = when the clusters are 4 Mpc apart, and 
central passage occurs at t ~ 2.15 Gyr. Approximately 0.5 Gyr 
after central passage a configuration is reached in which the overall 
gas morphology approximates fairly well the shape of the observed 
cluster. At this moment, the subcluster has passed beyond the major 
cluster centre while remaining significantly denser, which makes 
the emission slightly more intense in the region of the nose of the 
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http://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/ 
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Figure 1. Temporal evolution of model 233: (a) projected X-ray surface 
brightness, in the left column; and (b) projected emission-weighted tem- 
perature, in the right column. The six rows displays instants 0.2 Gyr apart. 
Each frame is 1.6x1.6 Mpc, the spatial scale is the same for all frames, and 
the surface brightness and temperature ranges are both constant. 
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shock than in the gas behind it. At earlier times no such distinction 
is noticeable, whereas at later times the global shape is excessively 
elongated and a detachment develops between the subcluster core 
and the gas left behind. 

In this particular model, the subcluster has 1/6 of the main 
cluster's total mass, but the gas in its core is more centrally concen- 
trated, being denser by a factor of 4. The subcluster is also colder, 
having an initial central temperature of ~ 1 keV, while the major 
cluster has 4 — 5 keV. The heated gas ahead of the cold subcluster is 
visible as a bow shock in the right column of Fig.[T]and it is further 
discussed in Section |4!4l 

At the instant of best match (t = 2.625 Gyr for model 233) 
the system evidently lacks spherical symmetry. Nevertheless we 
measure the spherically averaged density profile, centred on the 
point of highest total density. This is located at the major clus- 
ter's dark matter peak, around which the bulk of the total mass 
still lies (see Section l4~5b . We obtain the radius r^oo — 1-5 Mpc, 
at which the mean density has dropped to 200 times the critical 
density p c . This radius also encompasses the minor cluster's dark 
matter peak, as well as all the gaseous structures discussed here. 
The mass enclosed within r2oo is M200 — 4.0 x 1O 14 M0. The 
integrated X-ray luminosity (within r2oo and in the energy range 
0.1-2.4 keV) is Lx — 4.3 x 10 44 erg/s if the emission comes 
solely from bremsstrahlung, or Lx — 6.8 x 10 44 erg/s if metals 
and line emission are also taken into account. 



4.2 Exploration of parameter space 

In order to attempt to constrain some of the merger parameters, 
we explored numerous different initial condition configurations, in 
search of a 'best-fitting' model. Such a model would have to simul- 
taneously satisfy, to an approximate degree, the following criteria: 
the overall gas morphology should be reminiscent of the X-ray ob- 
servation of A3376; the temperature should be in the appropriate 
observed range; the total mass and total luminosity should fall in 
the same order of magnitude as those inferred from observations. A 
possibly important additional criterium regarding the distance be- 
tween the two brightest cluster galaxies is discussed in Section l4~5l 
We take the BCGs separation to be a proxy for the separation be- 
tween the two dark matter peaks. While model 233 may not nec- 
essarily be the one that strictly optimises each single criterium, it 
is the one that provided the most acceptable compromise among 
them. It is of course not possible to rule out the existence of al- 
ternative combinations of parameters that might result in similarly 
acceptable models. We explored physically motivated ranges of pa- 
rameters (albeit in a limited number of combinations due to the 
computational cost) and, at least for these ranges, certain configu- 
rations may me reliably excluded. 

As far as morphological comparison is concerned, we refrain 
from attempts to implement quantitative algorithms. Instead we 
rely on visual inspection, which may not be the most objective ap- 
proach and conclusions drawn from it ought to be regarded care- 
fully. Nevertheless, this is an approach not without its merits. As 
exemplified by the age-tested practice of morphological classifica- 
tion of galaxies, visual inspection tends to yield surprisingly reli- 
able results, which are often difficult to reproduce algorithmically. 
The type of comparisons we carry out here fall within this effort 
of making approximate judgements of morphology by eye, tak- 
ing into account both the overall appearance and various detailed 
aspects. While pixel-by-pixel subtraction (or some variation of it) 
could provide quantitative figures, it could hardly be expected to be 



sensitive to the same number of visual cues and at the same level 
of flexibility that evaluation by eye can afford. 

Here we present a systematic comparison of a sample of the 
models we explored. This is meant to allow a comparison of the X- 
ray morphology of the models listed in Table [TJ Taking model 233 
to be the standard, Fig.|2]displays variations around that model. The 
surface brightness scale and spatial scale are the same as in Fig.[T] 
Each row in Fig. [2] displays variations of one given parameter: (a) 
mass ratio; (b) impact parameter; (c) initial relative velocity; (d) 
relative central gas densities; and (e) inclination. For each of these 
properties, four variants are given, one of them being the fiducial 
model itself. This highlights the individual effect of varying each 
parameter separately. With the exception of row (e), every model is 
shown with the same inclination i — 40° to allow for a fair compar- 
ison. All snapshots are shown at the same instant t — 2.625 Gyr, 
with the exception of row (c), because models with different initial 
velocities have substantially different time-scales. 



4.2.1 Mass ratio 

Row (a) of Fig. [2] shows the outcome of cluster mergers having 
mass ratios of 1/2, 1/4, 1/6 and 1/8. It is clear that major merg- 
ers result in a quite distinct morphology. Models having Mb /Ma 
in the range of 1/6 - 1/8 are to be preferred. iNeistein & Dekell 
(2008), based on the extended Press-Schechter formalism, provide 
estimates of the rate of mergers having mass ratio above a certain 
value at a given redshift. A cluster of model 233's M200 at A3376's 
redshift will have undergone one merger per ~ 4.8 Gyr having 
M B /M a > 1/6. 



4.2.2 Impact parameter 

Typical impact para meters in galaxy cluster mergers are of the or- 
der of a few 100 kpc ( ISarazirj2002l;lRickeill998l;lRicker & Sarazinl 
l200lh . The four models shown in row (b) of Fig.|2]have initial im- 
pact parameters bo — 0, 150, 350 and 500 kpc. The most evident 
effect of off-centre mergers is the loss of symmetry around the col- 
lision axis. The curved trajectory of the subcluster is partly respon- 
sible for the distinctive shape of the resulting objects. In a head-on 
collision, the major cluster's central gas is spread out as the denser 
subcluster passes through it, decreasing its density. If however the 
two clusters don't pass through each other's centre, the major clus- 
ter's core remains relatively undisturbed, and as a result two sep- 
arate clumps of dense gas are still visible. The asymmetry due to 
a curved trajectory could be hidden from view if the orbital plane 
were seen exactly edge-on. But even under that particular configu- 
ration, both cores would still be discernible. Since A3 376 displays 
neither noticeable asymmetry nor two cores, there is no reason to 
go beyond the scenario of a head-on collision. 

It should be noted that bo refers to a shift in the direction of y- 
axis in the initial conditions, i.e. when the clusters are 4 Mpc apart 
in the a>axis direction. The minimum separation b m i n , which is 
the distance between the cluster centres at the instant of closest ap- 
proach, is considerably smaller. For the three off-axis models pre- 
sented in Fig.|2p, the approximate distances at pericentric passage 
are respectively b m i n — 100, 200 and 250 kpc. This sets a tight 
constraint, since a 6 mln as small as 100 kpc would be sufficient to 
produce noticeable asymmetry. 




Figure 2. Projected X-ray surface brightness for different models. Each row displays four variations of one given parameter: (a) mass ratio; (b) impact 
parameter; (c) initial relative velocity; (d) relative gas concentrations; (e) inclination. 
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Figure 3. Comparison between observations of A3376 (left) and model 233 
(right). The upper row displays the XMM-Newton X-ray surface brightness 
in the [0.5-8.0 keV] band compared to a simulated image whose resolu- 
tion has been deliberately degraded. In the lower row, the left panel shows 
the observed temperature maps. In the simulated temperature map, a semi 
transparent mask highlights the region in which there is data. 

4.2.3 Initial relative velocity 

Given the major cluster mass Ma, the subcluster's free fall velocity 
at do = 4 Mpc would be approximately vo = 1250 km/s, if it were 
a point mass having been released from infinity at rest. 

The models displayed in row (c) of Fig. [2] have initial rela- 
tive velocities vo = 500, 1000, 1500 and 2000 km/s. Because the 
time-scales are not the same, they are compared at different in- 
stants, each chosen to be that which best resembles A3376. As far 
as the gas morphology is concerned, the different velocities affect 
the sharpness of the nose of the shock. Models vlOOO and v2000 
were dismissed not on account of their X-ray morphologies, which 
are not inadequate. However, model v2000 gives rise to excessively 
high temperatures in the shock region. In model vlOOO, on the other 
hand, it takes a long time for the dark matter peaks to be sufficiently 
separated (see Section |4~5t and by that time the morphology has de- 
teriorated. The 1500 km/s model provides a tolerable compromise 
between morphology, temperature and dark matter peaks separa- 
tion. The Mach numbers corresponding to these velocities are dis- 
cussed in Section l4"l4l 

4.2.4 Central gas density 

Row (d) of Fig.|2]highlights the effects of relative central gas densi- 
ties, i.e. the effects of the concentration of the subcluster. It shows 
models with no t B/no,A = 1, 2,4 and 6. There is a clear depen- 
dence on the concentration of the subcluster. When the clusters 
have comparable concentrations, the subcluster gas is hardly distin- 
guishable from the major cluster gas. The outcome is a somewhat 
uniform emission with no X-ray peak. If the subcluster is consider- 
ably denser, it is able to cross the major cluster core and to remain 
relatively cohesive as it emerges. Because in model n6 the X-ray 
peak is excessively prominent, the preferred model is that in which 
the subcluster central gas is 4 times denser the major cluster's. 

4.2.5 Inclination 

Row (e) of Fig. [2] shows model 233 in the same instant projected 
under four different inclinations i — 0°, 30°, 40° and 60°. For low 
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Figure 4. The upper panel displays the shock velocity v s (dot-dashed line) 
and the velocity of the upstream gas u (dashed line), both measured in the 
centre-of-mass rest frame; also shown is the sound speed c s (solid line). 
The lower panel shows the resulting Mach number as a function of time for 
this model (233). The shaded regions represent upper and lower boundaries 
for each quantity, computed from models vlOOO and v2000 (and shifted to 
match the same time interval). 



inclinations, the shape is excessively elongated at this time. For 
very high inclinations it is excessively round. Unfortunately the 
observations provide no information on the spatial orientation of 
the cluster. Here again, the best configuration is chosen not solely 
from the X-ray morphology, although it is a good indicator. For 
example, at later time an inclination of i = 60° provides a mor- 
phology that is also acceptable. However the projected separation 
of the dark matter peaks is very time dependent and also very sen- 
sitive to inclination (see Section |4~5l l. By the time the dark matter 
peaks are sufficiently separated with i = 60°, the morphology is 
no longer adequate. Therefore the i — 40° projection provides the 
best compromise. 



4.3 Comparison to observations 

A comparison between observations of A3376 and model 233 is 
given in Fig. [3] Despite the high-resolution of the numerical model 
itself, the simulated images are deliberately degraded by undersam- 
pling the particles, by applying a gaussian smoothing and by the 
addition of noise. 

In the observed temperature map (see Section 13.21 there is 
data only within a relatively small region approximately 20 ar- 
cmin wide. To make the comparison more straightforward, a re- 
gion of the same shape is overlaid on the simulated temperature 
map, while the surrounding area is made semi opaque. This empha- 
sises the region in which data exists, and underscores the simulated 
temperature features which are not observationally available. Apart 
from the range of values, the observed temperature map does not 
exhibit remarkable features that could impose particularly strong 
constraints on the simulations parameters. We were able to rule 
out collisions that heated the gas to above 10 keV, for example, 
and this was used to constrain velocities and concentrations. Nev- 
ertheless, the small scale details of the observed temperature map 
are not reproduced, probably due to the simplifying assumptions 
of our adiabatic numerical models, which take into account neither 
substructures nor the various galaxy-related physical processes that 
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might affect the intracluster gas, such as feedback from supernovae 
and AGN. 



4.4 Mach number 

In a cluster merger simulation, the Mach number can be directly 
measured, as all velocity information is available. One of the most 
pronounced features is the temperature drop after the bow shock. 
The successive positions of this discontinuity are used to directly 
compute the shock velocity v s ( in the centre-of-mass rest frame). 
As shown bv lSpringel & Farrarl (12007b in simulations of 1E0657- 
56, the gas ahead of the shock is not at rest. In such mergers, the 
upstream gas is in fact falling towards the incoming subcluster with 
a considerable velocity u (in the centre-of-mass rest frame). There- 
fore, the effective relative velocity with which the shock front en- 
counters the pre-shock gas is v 3 — u. Consequently, the Mach num- 
ber ought to be computed as 



M 



(5) 



where c s — is the sound speed of the pre-shock gas. Both the 
upstream velocity and the sound speed are measured in the region 
immediately ahead of the shock. 

For model 233 at t = 2.625 Gyr the shock front moves for- 
ward with v s = 2114 km/s while the upstream gas falls back 
at u = —468 km/s. The sound speed c s = 666 km/s implies 
A4 = 3.9 at this instant. Figure [4] shows how these quantities 
evolve over an interval of 0.4 Gyr. To provide an indication of the 
typical ranges of these velocities, the shaded areas represent the 
boundaries given by models vlOOO and v2000 (shifted to match 
their respective time-scales to that of model 233). The resulting 
Mach numbers, bounded by these extreme cases, would be in the 
range of roughly M — 3 — 4. 

Once A4 has been computed directly from the velocities, it 
may be used to obtain the expected shock discontinuities from 
the Rankine-Hugoniot conditions. Figure[5]displays five quantities 
measured along the collision axis: temperature; electron number 
density; pressure; entropy; and gas streaming velocity in the x- 
axis direction. As a proxy for entropy, the conventional definition 
S = k T nj 2 ^ 3 is adopted. All of these profiles have been mea- 
sured using not the projected images, but using the particles within 
a cylinder of radius 150 kpc passing through the nose of the shock. 
The upper panels of Fig. [5] display the surface brightness and the 
temperature both projected under i — 0° inclination. The vertical 
dashed lines at x — 1148 kpc mark the position of the shock front, 
which has been determined as the point where the temperature drop 
is most intense. The vertical solid lines at x — 980 kpc mark the po- 
sition of the contact discontinuity (the cold front), the point where 
the density drops the most. At the contact discontinuity velocity and 
thermal pressure are continuous. The horizontal lines indicate the 
expected drops in each of the five quantities at the shock position, 
as computed from the Rankine-Hugoniot conditions using the mea- 
sured A4 . Assuming 7 = 5/3 throughout, the relations between the 
pre- shock (subscript 1) and post-shock (subscript 2) quantities are 
(e.g. lLandau & Lif shit J 19591 ; [shiill 1 992h : 
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Figure 5. Properties of the shock for model 233 at t = 3.3 Gyr. The first 
and second panels show respectively the X-ray surface brightness map and 
the emission-weighted temperature map. The third to seventh panels display 
the following quantities measured along the a>axis: temperature, electron 
number density, pressure, entropy, and gas streaming velocity. The location 
of the contact discontinuity is marked by solid vertical lines; the location 
of the bow shock is marked by dashed vertical lines. The dashed horizontal 
lines give the values of each quantity before and after the shock, where 
the expected drops were computed from the Rankine-Hugoniot conditions 
using the measured Mach number. 
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Figure 6. Temperature and density measured along the x-axis for model 
233 seen in projection with inclination i = 40°. If the temperature drop 
is measured from the inclined model, the inferred Mach number is smaller 
than the actual value. 
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The good match between the expected drops and the measured 
profiles indicates that the Mach number could, in principle, be in- 
ferred from these quantities. For an observed shock, there is how- 
ever the difficulty introduced by the unknown inclination. To illus- 
trate this problem, we take the same model 233 at the same instant 
in time and now project it under the inclination i — 40°. If we 
now try to measure the temperature and density along the projected 
collision axis, we obtain the profiles shown in Fig. [6] The result is 
that the height of the temperature peak is lowered because, under 
projection, the thin region of very hot gas is seen as spread over 
a larger area. Furthermore, the sharpness of the discontinuities is 
attenuated. Now, instead of using the known Mach number to com- 
pute the drops, we conversely measure the temperature and density 
drops (horizontal lines in Fig. [6] at the shock position. These ratios 
would lead to A4 = 2.9 for this shock. A Mach number inferred 
in this manner for a cluster of unknown inclination should thus be 
regarded as a lower limit. 



4.5 Dark matter distribution 

In A3376 the brightest cluster galaxy (first BCG) is located not 
at the peak of X-ray emission, but at a distance of approximately 
970/ifg 1 kpc from it. The X-ray peak coincides with the second 
BCG. A possible interpretation of this peculiarity is a scenario in 
which a less massive cluster comes from the southwest hosting the 
second BCG and overtakes the major cluster's core, where the first 
BCG lies. Our best model accounts for this feature, in the following 
sense. As there are no galaxies in our simulations, we assume the 
BCG separation may be identified with the dark matter peaks sepa- 
ration. In the absence of additional clues, it is reasonable to assume 
that the most massive galaxies should in principle be located at the 
bottoms of the potential well, i.e. at or near the centroids of the two 
dark matter peaks. 



Figure 7. The upper panel shows the DSS optical image in which the circles 
mark positions of the galaxies at the cluster redshift. The first BCG (right) 
and the second BCG (left) are highlighted. In the lower panel, white coun- 
tour lines show the total projected mass, overlaid on the projected X-ray 
surface brightness, for model 233 at t = 2.625 Gyr with i = 40°. 



The lower panel of Fig.|7]shows white contour lines that rep- 
resent the total projected mass for model 233, overlaid on the X-ray 
surface brightness. The separation between the dark matter peaks 
is approximately 800-850 kpc at the best-fitting instant. The upper 
panel of Fig. [7] shows an optical 'true-colour' (IR-Red-Blue) DSS 
image of A3376 in which the galaxies at the cluster redshift (se- 
lected using NEE0) are marked by circles and the two BCGs are 
highlighted. The dark matter separation of model 233 matches the 
observed BCG distance to within ~ 12%. Again, an even better 
match would be achieved at a slightly later time, but at the price of 
a deteriorating gas morphology. 

Of course, a possible mass map of A3376 exhibiting two ma- 
jor mass concentrations around the locations of the BCGs would 
lend more credence to this scenario. Mapping the dark matter distri- 
bution by means of gravitational weak lensing is particularly chal- 



NASA/IPAC Extragalactic Database, http://ned.ipac.caltech.edu/ 
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lenging in the case of A3376 due to its proximity. Once or if such 
map is available, it might either corroborate this scenario or over- 
throw it. If the dark matter separation turns out to be similar to the 
BCG separation, the mass map might set more stringent constraints 
on the simulation parameters. If, on the other hand, a more compli- 
cated dark matter structure is revealed, then alternative models will 
have to be sought out. 

An additional degree of freedom that has not been explored 
in depth in this paper is the relative concentration of the dark mat- 
ter haloes. For example, if the dark matter is too centrally concen- 
trated, gas temperatures need to be in excess of 10 keV even in the 
initial conditions to satisfy hydrostatic equilibrium. To obtain phys- 
ically plausible temperatures of 5 keV in the major cluster, we em- 
ploy halo scale lengths Th ~ 500 kpc comparable to the gas density 
scale lengths r g . If the dark matter concentrations are too low, the 
subcluster's dark matter is able to advance further than its gas and 
a dissociation develops. As of yet, there is no reason to believe that 
this is the case for A3376. If a dark matter/gas offset is shown to 
exist, then it would have to be accommodated by exploring differ- 
ent dark matter concentrations. A more systematic analysis of the 
dark matter distribution and its dependence on merger parameters 
is beyond the scope of this paper. 



5 SUMMARY AND CONCLUSIONS 

The peculiar cometary shape of A3376 is suggestive that it has 
undergone a recent merging event. This is, to best of our knowl- 
edge, the closest galaxy cluster displaying such morphology. Fur- 
thermore, the diffuse radio emission in the form of double Mpc- 
scale arcs in its periphery is un derstood to be cau sed by shock- 
accelerated relativistic electrons (Bagchi etal. 200 j). 

Using cosmological simulations, iPaul et al. obtains a 

cluster whose shock wave structure resembles the radio relics of 
A3376. Dedicated hydrodynamical simulations meant to model one 
specific merging clus ter have been mostly focused on the Bul- 
let Cluster itself (e.g. iMilosavlievic et al.ll2007l : ISpringel & Farrarl 
I2007L iMastropietro & Burked 120081) . Recently Ivan Weeren et al.l 
( 201 if) and Briiggen et al. (2012) carried out simulations of the 
galaxy clusters CIZA J2242.8+5301 and 1RXS J0603.3+4214, re- 
spectively, both of which also exhibit radio relics. In the case of 
A3376, the morphology is sufficiently uncomplicated that it may 
be satisfactorily modelled as the encounter of only two objects, thus 
rendering the reconstruction of its dynamical history relatively sim- 
pler. 

The simple numerical model we set up in order to investigate 
the problem consists in the collision of two spherically symmet- 
ric galaxy clusters. Equilibrium initial conditions are prepared in 
which the galaxy clusters are represented by a gas component and 
a dark matter component. The hydrodynamical simulations them- 
selves are adiabatic, and we assume that radiative losses are unim- 
portant during the time span of the simulations. Yet we use the out- 
put of the simulations to generate maps of projected X-ray surface 
brightness and maps of emission-weighted temperature, in order to 
compare them to XMM data. 

Starting from a separation of 4 Mpc and a relative initial ve- 
locity of 1500 km/s, the clusters meet on a head-on collision and 
the best-matching moment is reached approximately 0.5 Gyr af- 
ter central passage. In order to constrain some of the collision pa- 
rameters, we covered as much parameter space as allowed by the 
computationally intensive nature of such effort. Ideally, in order 
to reach a 'best model' a number of criteria would have to be si- 



multaneously met, namely: overall gas morphology, temperature, 
virial mass, total X-ray luminosity, and distance between the two 
dark matter peaks. While the best model presented here may not 
necessarily optimise each of these criteria individually, it provided 
the most adequate compromise. It is of course impossible to argue 
for the uniqueness of a solution, as alternative combinations of pa- 
rameters could conceivably provide similar outcomes. This impos- 
sibility notwithstanding, the physically motivated ranges of values 
explored here at least allow us to reliably rule out certain combina- 
tions of parameters. 

Here we summarise the five main parameters that have been 
constrained and give approximate estimates of the best ranges: (a) 
We find that the best matches are obtained for mass ratios in the 
interval 1/6 - 1/8. Major mergers are excluded due to their global 
morphology, while in minor mergers of very small mass ratio, the 
X-ray peak is not sufficiently intense, (b) Large impact parameters 
are straightforwardly ruled out as they lead to obvious asymme- 
tries that are not observed in A3376. An approximate upper limit 
to 60 is set at 150 kpc, which implies a separation of < 100 kpc at 
pericentric passage, (c) The collision velocity is constrained not by 
the morphology alone, but also by the temperatures and dark mat- 
ter peak separation. For initial velocities of 2000 km/s the resulting 
temperatures are excessively high in the shock region, whereas for 
1000 km/s the desired dark matter peak separation takes too long 
to be reached, (d) An important parameter is the relative central gas 
concentration. We find the most adequate morphology is obtained 
when the subcluster is denser than the major cluster by a factor of 
about 4 in the centre. This determines essentially the prominence of 
the X-ray peak, which is excessively outstanding if the subcluster 
central density is too high, or nearly unnoticeable otherwise, (e) Fi- 
nally, the inclination (the angle between the collision axis and the 
plane of the sky) plays an important role in the morphology and 
it is strongly time-dependent. Both the temperature and the total 
X-ray luminosity are greatly increased at central passage. By the 
time they reach adequate levels, the shape of the gas distribution 
is excessively elongated when viewed on the collision plane. We 
find that, projected under an inclination angle of i = 40°, the mor- 
phology is considerably less elongated and the separation between 
the two dark matter peaks is comparable to the BCG separation to 
within ~ 12%. 

Merging clusters generally exhibit complicated temperature 
structures. The observed temperature map, in a region approxi- 
mately 20 arcmin wide, shows a mean temperature of ~ 3.5 keV, 
but no discernible features that could be used to set strong con- 
straints on the simulations. We use it to rule out models in which 
the shock heats the gas considerably above the observed range. The 
outcome of the simulations suggests that the region in which there 
is data encompasses the contact discontinuity at most, but excludes 
the shock front itself, i.e. the shock-heated gas ahead of the sub- 
cluster's X-ray peak. Furthermore, the small scale details of the 
observed temperature map are not reproduced in our simulations. 
If they are the result of interactions between galaxies and the intr- 
acluster medium, they could not have been recovered by our sim- 
plified simulations which include no such physical processes. The 
simplifying assumptions of these simulations also exclude the ef- 
fects of additional substructure. 

Cluster mergers are expected to d rive supersonic shock waves 
of typical Mach numbers M < 3 dSarazinl 120021) but stronger 
shocks may arise under som e circumstances (e.g. IVazza et all 



1201 ll : IPlanelles & Ouilisll2012l). From Suzaku X-ray observations 
of A3376, lAkamatsu et all d2012h were able to measure a tempera- 
ture jump in the western radio relic leading to A4 = 2.91 ± 0.91. 
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In our simulations, a supersonic shock wave develops ahead of the 
colder subcluster. The peak of X-ray emission corresponds to a low 
temperature, dense gas ahead of which the shock front is located. 
The layer of heated gas between the edge of the subcluster and the 
bow shock is ~ 170 kpc thic k. We determine the shock velocity 
as in[Springel & Farrar (2007), taking into account the velocity of 
the upstream gas that is falling towards the incoming subcluster. 
The actual shock velocity, i.e. the velocity with which the shock 
wave meets the pre-shock gas, is ~ 2600 km/s, resulting in a mach 
number of 3.9. This was computed from the intrinsic properties of 
the simulation output. If, alternatively, we use the projected im- 
ages (with inclination i = 40°) to measure the apparent jumps in 
the temperature and density profiles, we obtain M — 2.9 which 

should then be regarded as a lower limit. 

Besides the notorious bullet cluster dClowe etal ] |2006l) , there 
is a growing list of merging clusters that have been shown 
to have an offset between their dark matter and their gas 



(e.g. iMahdavi et alj|2007l; iBradac et alj|2008l ; lOkabe etalj|201ll ; 
Daws onet all2012l ; lRagozzine et alj|2012l) . In our best fitting sim- 
ulations, no substantial dissociation between gas and dark matter 
developed, and no such models were pursued because, as of yet, 
there is no observational evidence that this type of offset took place 
in A3376. Since the total mass distribution of A3376 is unknown, 
we use the locations of the two brightest cluster galaxies as indi- 
cators of the positions of the dark matter peaks. In the absence of 
further data, it is reasonable to assume that the BCGs should in 
principle coincide with the centroids of the two dark matter haloes. 
In this scenario, the major cluster is assumed to host the first BCG, 
and the subcluster is believed to carry the second BCG, whose posi- 
tion coincides with the peak of X-ray emission. If this is so, the bulk 
of the dark matter in the system should be located ~ 970/if 1 kpc 
away from the X-ray peak. This is what we find in our simulations: 
a secondary dark matter peak at the position of the X-ray peak, and 
the main dark matter peak ~ 850 kpc behind it. The dynamics of 
how the two dark matter haloes go through each other depends on 
their relative concentrations and has relevant effects on the final lo- 
cations of the dark matter peaks. However, a systematic analysis of 
this aspect will be presented elsewhere. 

Determining the projected mass of A3376 by the technique of 
gravitational weak lensing is difficult because the cluster distance 
is relatively small. If this proves to be feasible, the resulting mass 
map could conceivably set tighter or additional constraints on the 
current merger parameters. If an unexpected dark matter distribu- 
tion is uncovered, then the current scenario might turn out to be 
insufficient and a further exploration of the least robust parameters 
will be required. 

Bearing in mind the limitations of this approach by TV-body 
simulations, we have proposed a specific scenario of a dynamical 
history for the merging event of A3376 and offered a possible com- 
bination of parameters that accounts for several of its features. Fu- 
ture weak lensing analysis might either help corroborate this picture 
or necessitate its improvement. 
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